Fire Severity Analysis

Fire Severity Analysis in Guatemala Using Landsat 8/9 Imagery

1 Introduction

Tutorial for calculating and plotting dNBR Fire Severity Index using Landsat 8/9 satellite imagery.

1.1 Data Sources

USGS Earth Explorer

1.2 Sources Referenced

The code in this analysis was adapted from the following sources below:

Wildfire Severity Analysis

Landsat Remote Sensing tif Files in R

1.3 Additional Contributors

Conservation Data Lab, Randy Swaty, Sam Ettenborough

2 Setup

Please review the README.md document for installation and setup instructions.

3 Importing Landsat Files

# GEE project title
project <- "ee-samettenborough"

# fire 1 coordinates
fire_1 <- list(-91.253, 16.8194, -89.4018, 18.2112)

# these two fires started end of may and last measurements were in end of july
aoi_03 <- here::here("data/areas-of-interest/EMSR727_AOI03_BLP_PRODUCT_areaOfInterestA_v1.shp")
aoi_05 <- here::here("data/areas-of-interest/EMSR727_AOI05_BLP_PRODUCT_areaOfInterestA_v1.shp")

# pre fire start and end date search window
pre_start_date <- "2024-01-01"
pre_end_date <- "2024-03-01"

# output folder
out_dir <- "data/pre"

# uncomment this section to re-download the images
# process_aoi(
#     project,
#     pre_start_date,
#     pre_end_date,
#     out_dir,
#     geometry = aoi_05,
#     prefix="aoi_05_landsat_pre_SR5_SR7_",
# )


# pre fire start and end date search window
post_start_date <- "2024-04-01"
post_end_date <- "2024-08-01"

# output folder
out_dir <- "data/post"

# # uncomment this section to re-download the images
# process_aoi(
#     project,
#     post_start_date,
#     post_end_date,
#     out_dir,
#     geometry = aoi_05,
#     filter_lim = 30,
#     prefix="aoi_05_landsat_post_SR5_SR7_",
# )

4 Visualize Pre and Post Images

pre_fire_files <- list.files("data/pre/",
                             pattern = "aoi_05_landsat_pre_SR.*",
                             full.names = TRUE)

pre_fires <- pre_fire_files |>
  map(rast) |>
  reduce(mosaic)

# Plotting in RGB to test
par(col.axis = 'white', col.lab = 'white', tck = 0)
plotRGB(pre_fires,
        r= 2, g= 1, b= 1,
        stretch = 'lin',
        axes = TRUE,
        main = 'Prefire RGB Composite Image, Bands 5,7')

# # Create leaflet map
# leaflet() %>%
#   addTiles() %>%
#   addRasterImage(
#     pre_fires, 
#     opacity = 0.7,
#     group = "Pre-Fire RGB Composite"
#   )

post_fire_files <- list.files("data/post",
                              pattern = "aoi_05_landsat_post_SR.*",
                              full.names = TRUE)

post_fires <- post_fire_files |>
  map(rast) |>
  reduce(mosaic)

par(col.axis = 'white', col.lab = 'white', tck = 0)
plotRGB(post_fires,
        r= 2, g= 1, b= 1,
        stretch = 'lin',
        axes = TRUE,
        main = 'Postfire RGB Composite Image, Bands 5,7')

5 Calculate Normalized Burn Ratio and Burn Severity (dNBR)

# Calculating the Delta Normalized Burn Ratio Index 
nbr_pre <- 1000 * (pre_fires[[1]] - pre_fires[[2]]) / 
  (pre_fires[[1]] + pre_fires[[2]])

nbr_post <- 1000 * (post_fires[[1]] - post_fires[[2]]) / 
  (post_fires[[1]] + post_fires[[2]])

# Attempt to calculate the Differenced Normalized Burn Ratio Index
dnbr <- (nbr_pre)-(nbr_post)

6 Display Pre and Post Fires

nbr_stack <- c(nbr_pre, nbr_post)
names(nbr_stack) <- c("Pre-fire NBR", "Post-fire NBR")
nbr_stack_df <- rasterdf(nbr_stack)

ggplot(nbr_stack_df) +
  geom_raster(aes(x = x, 
                  y = y, 
                  fill = value)) + 
  scale_fill_gradient(name = "NBR", 
                       low = "lightyellow", 
                       high = "darkgreen") +
  coord_sf(expand = FALSE) +
  annotation_scale(location = 'bl') +
  facet_wrap(facets = vars(variable), 
             ncol = 1) +
  theme_void()

7 Display dNBR

dnbr_df <- rasterdf(dnbr)

ggplot(dnbr_df) +
  geom_raster(aes(x = x, 
                  y = y, 
                  fill = value)) + 
  scale_fill_gradient2(name = "DNBR", 
                       low = "blue", 
                       high = "red",
                       midpoint = 0) +
  coord_sf(expand = F) +
  annotation_scale(location = 'bl') +
  theme_void()

8 Classify Severity

# Classifying the index values into a matrix
rclas <- matrix(c(-Inf, -970, NA,  # Missing data
                  -970, -100, 5,   # Increased greenness
                  -100, 80, 1,     # Unburned
                  80, 265, 2,      # Low severity
                  265, 490, 3,     # Moderate severity
                  490, Inf, 4),    # High severity
                ncol = 3, byrow = T)

severity <- classify(dnbr, rclas)

SCcolors = c("#008080", 
             "#5f9ea0", 
             "#e0e0e0", 
             "#a0522d", 
             "#8b0000")
SCnames = c("Unburned", 
            "Low", 
            "Moderate", 
            "High", 
            "> Green")

severity_df <- rasterdf(severity)

ggplot(severity_df) +
  geom_raster(aes(x = x, 
                  y = y, 
                  fill = as.character(value))) + 
  scale_fill_manual(name = "Severity Class",
                    values = SCcolors,
                    labels = SCnames,
                    na.translate = FALSE) +
  annotation_scale(location = 'bl') +
  coord_fixed(expand = F) +
  theme_void()

severity_counts <- severity_df %>%
  count(value) %>%
  mutate(
    percentage = (n / sum(n)) * 100,  # Convert counts to percentages
    value = as.character(value)       # Ensure value is a character for plotting
  )

# Create the bar chart showing percentages
ggplot(severity_counts, aes(x = value, y = percentage, fill = value)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(name = "Severity Class", values = SCcolors, labels = SCnames) +
  labs(x = "Severity Class",
       y = "Percentage of Pixels",
       title = "dNBR Percentages by Severity Class",
       subtitle = "A0103: Paso Caballos") +
  theme_minimal() +
  geom_text(aes(label = sprintf("%.1f%%", percentage)),  # Format label to 1 decimal place
            vjust = 1, size = 5) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),  # Center title & increase size
    axis.title = element_text(size = 14)
  )

severity_counts <- severity_df %>%
  count(value) %>%
  mutate(
    percentage = (n / sum(n)) * 100,  # Convert counts to percentages
    value = as.character(value)       # Ensure value is a character for plotting
  ) |>
  filter(!is.na(value))

# Create the bar chart showing percentages
severity_plot <- ggplot(severity_counts, aes(x = value, y = percentage, fill = value)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(name = "Severity Class", values = SCcolors, labels = SCnames) +
  labs(x = "Severity Class",
       y = "Percentage of Pixels",
       title = "dNBR Percentages by Severity Class",
       subtitle = "A0105: Sierra del Lacandon") +
  theme_minimal() +
  geom_text(aes(label = sprintf("%.1f%%", percentage)),  # Format label to 1 decimal place
            vjust = 1, size = 5) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),  # Center title & increase size
    plot.subtitle = element_text(hjust = 0.5, size = 14)
  )

ggsave("outputs/A0105-dNBR-Percentages-by-Severity-Class.jpg",
       plot = severity_plot, width = 8, height = 6, dpi = 300)

# # Define severity classification matrix
# rclas <- matrix(c(-Inf, -970, NA,   # Missing data
#                   -970, -100, 5,     # Increased greenness
#                   -100, 80, 1,       # Unburned
#                   80, 265, 2,        # Low severity
#                   265, 490, 3,       # Moderate severity
#                   490, Inf, 4),      # High severity
#                 ncol = 3, byrow = TRUE)

# Color palette
severity_colors = c("#008080", 
             "#5f9ea0", 
             "#e0e0e0", 
             "#a0522d", 
             "#8b0000")
severity_names = c("Unburned", 
            "Low", 
            "Moderate", 
            "High", 
            "> Green")

# Create leaflet map
leaflet() %>%
  addTiles() %>%
  addRasterImage(
    severity,
    colors = severity_colors,
    opacity = 0.7,
    group = "Severity Classes"
  ) %>%
  addLegend(
    position = "bottomright",
    colors = severity_colors,
    labels = severity_names,
    title = "Severity Class"
  )